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Abstract 

We apply belief propagation (BP) to multi-user detection in a spread spectrum 
system, under the assumption of Gaussian symbols. We prove that BP is both 
convergent and allows to estimate the correct conditional expectation of the input 
symbols. It is therefore an optimal -minimum mean square error- detection algo- 
rithm. This suggests the possibility of designing BP detection algorithms for more 
general systems. 

As a byproduct we rederive the Tse-Hanly formula for minimum mean square 
error without any recourse to random matrix theory. 

1 Introduction 

Consider the multiuser detection problem of K users, each spreading its symbol 
onto N chips in a spread spectrum system 

K 

i=i 

Here Xi € R is the symbol transmitted by user i, Si G M> N is the iV-chip long 
signature sequence of user i, y is the received vector, and w is an N x 1 vector 
of i.i.d. mean and variance a 2 Gaussians. We assume that the Xi are i.i.d. with 
Kxi = and Ei ( 2 = 1. Finally, we shall denote by § the N x K signature matrix 
whose columns are the vectors s,. 

Several years ago, Tse and Hanly 1 , and Verdu and Shamai [2] considered the 
case in which the symbols Xj's are Gaussian random variables, which models the 
situation when they are ideal Shannon-coded symbols. The signature sequences 
were assumed themselves random, but known both at the receiver. More precisely 
(here and in the following denotes the transpose of matrix, or vector, ^4) 



Si = —=(sn, s iN ) T , for i = 1, K, 
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Figure 1: Left: A graphical representation of the multiuser communication problem. 
Right: Messages notation in the BP detection algorithm. 

where the Sj a 's are i.i.d. random variables with zero mean and unit variance. These 
authors considered a minimum mean square error (MMSE) single user receiver 
which generates soft estimates Xj's and analyzed its performance in the large system 
limit N,K — > oo, with a = K/N constant. The analysis was heavily based on 
random matrix theory and, in particular, on the Marcenko-Pastur theorem on 
the eigenvalue distribution of large random matrices. As a consequence, it did 
not suggest any generalization to the practically interesting case of non-Gaussian 
symbols Xj's. 

Recently, Tanaka [5j applied the replica method from statistical mechanics to 
the case of uncoded binary antipodal signals: xi S {+1,-1}. He was able to 
compute the asymptotic bit error rate and conditional entropy per bit, in the large 
system limit. Given the relationship between the replica and cavity method from 
statistical mechanics and message passing algorithms, such as belief propagation 
(BP) and various generalizations thereof, several authors [U [3] studied the use of 
BP as a detection algorithm in this context. If proved to be correct, this approach 
would provide a low complexity, asymptotically optimal receiver. 

Although the results obtained via replica method are likely to be correct 1 , the 
method itself is non-rigorous. Also, the studies of BP detection algorithm were 
essentially empirical. No convergence or correctness guarantees exist. 

In this paper we follow the opposite route to the one sketched above. We 
consider BP detection and rigorously prove convergence and correctness. In a 
second step, we analyze BP to compute the system performance. This approach 
solves both the problems stressed above (non rigorous character of the replica 
method, and lack of guarantees on BP detection). This paper deals with the case 
of Gaussian symbols x,'s, with zero mean and unit variance. Following our strategy, 
we are able to recover several of Tse and Hanly's results without any recourse to 
random matrix theory. We conjecture however that the same approach can be 
developed for non-Gaussians symbols allowing to recover Tanaka's results without 
any recourse to the replica method. 

In order to apply BP, it is convenient to formulate multiuser detection as an 
inference problem on a probabilistic graphical model. The underlying graph is 
depicted in Fig. ^ It is a complete bipartite graph on K left nodes (users) and N 
right nodes (chips). We associate real variables x = (x\, . . . , xk) to the user nodes 

1 Tanaka's results are derived under the 'replica symmetry' assumption. There are arguments sup- 
porting this assumption in the problem as described so far. On the other hand, the same hypothesis 
needs to be replaced by more refined ones for generalizations of this problem (for instance if the noise 
variance is not known at the receiver). 



and u = (o>x, • • • ,^n) to the chips, and consider the (complex) weight 
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d^ K (x,ui) = ^ m ll e-3^+^J]di/(x i )nexp|--^ So ia; a x i | dw. (3) 
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where d^(xj) = exp(— ^xf)/2ir is the a priori distribution of the symbol Xj, and 
j = \/— T. Elementary calculus shows that dfj,y' K (x) = f dfj,y' K (x,u) is in fact 
the conditional distribution of the transmitted symbols x given y. The detection 
problem amounts to compute the marginal dfj, y ' (xi), i.e. to integrate Eq. (J3J) 
over all the variables but X{. Since, the weight (jSJ) is Gaussian, d/j,y ,K (xi) will be 
Gaussian as well and can be parameterized as 

dfiy' K (xi) = exp j-^Lixf + GixX dx , (4) 

The weight @ factorizes according to the complete bipartite graph in Fig. ^ left 
frame. BP can therefore be used to computing the marginals d^Ly ,K ' (x{) through a 
message passing procedure. Messages are exchanged as in Fig. Q right frame, and 
computed according to the update equations 

x(*+l) I I 1 s ib TW _ _2 , 1 S la ( r\ 

with initial conditions A-^ a = 1, 7 4 -__^ a = 0. The parameters in the marginal 
distribution d[i y ' K (xi), cf. Eq. @, are estimated as 

n(t+l) _ 1 S ib fJt) j r(t+l)_-| , 1 4 m 

G * -^Z^TOT^, and ^ _i+ ^^r- (7) 
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In the next Section we will show that this procedure is convergent and provides 
the asymptotically correct marginal distributions, thus implementing allowing to 
implement a MMSE estimator. In Section |3] we show that the usual prescription 
for a MMSE receiver can be recovered from the fixed point of BP. Finally, in 
Section^ we investigate the rate of convergence of our algorithm through numerical 
simulations. 



2 Main result: convergence and correctness 

Throughout this Section we assume that Si a = ±1 with equal probability. We 
believe that our main results, as well as the idea of the proof, remain valid for a 
considerably more general distribution. However, this assumption allows to avoid 
several technical complications. In particular, it is immediate to prove the following 
result 

Lemma 1 If s ia G {+1,-1} for any i G [K], a G [N], then \fi a = X(t), X^Li = 
X(t), where X(t),X(t) are given by the iterations 



(8) 



with A(0) = 1. 

Moreover X(t) -> \(°°\ X(t) -> X^ as t -> oo, with A(°°\ A^) > the unique 
positive fixed point of the above equations. More precisely, we have A^°°^ = 1 + A, 
A(°°) = 1/A, where 

1 9 a 

a=" + TTa' < 9 > 

The reader will recognize that Eq. © is nothing but the Tse-Hanly asymptotic 
formula for the signal to interference ratio in the large system limit. 

We shall state two separate results for the mean of the marginal distribution 
Q, Xi = Gil Li and its variance 1/Li. 

Theorem 1 Assume the Si a 's are i.i.d. uniformly random in {+1,-1}. Then, 
there exists a set S N C {+l,-l} KxN of signatures with F(S N ) > 1 - 0(N~ a ) 
(with a an arbitrary positive number) such that, for any 8 £ Sn and for every 
{y a }, L^ — > a^c? Gf — > g[°°^ as t — > oo. Furthermore, the conditional 

means estimated by Belief Propagation are correct, i.e., x 

(oo) = G (oo) /L (oo) 

Sketch of proof: For Gaussian graphical models, it was proved in [3 |H] that, 
if BP converges for a 'generic initial condition', then it correctly estimates the 
mean of the associated probability distribution. Referring to the present case, this 
immediately implies = Xi if convergence holds for any initial condition of the 

form 7 £2 a = Ci^ a , Xf_l a = 1- 

Because of Lemma^ it is in fact sufficient to prove convergence for the messages 
{7j-^. a }- We consider therefore the homogeneous recursion 

We shall show that for a typical realization of the signatures S, and rjj*^ — > 
for any initial condition r|^ a and r a °^. This in turns imply that ^if\, a and 
converge. We begin by eliminating from the above relations to get 



vfzj = 7T-7T E n ^ b r^ b (ii) 



where 



JO if i = A; or a = b, . . 

[ •SjfcSfcb otherwise. 

The random matrix f2 = {Oj a ^ : i,k E [X], a,b € [A]} has dimensions Af-fiT x A^i? 
and is a function of the signature sequence S. 

We make use of the following property of $7, whose proof is omitted. 

Lemma 2 For each a > 0, and each k > 1 there exist a positive integer Nq = 
N (a,n) such that, if N > N and t < (AylOOO) 1 / 6 then 

2t+2„,i+l ( i , C(a, k) 4 q 



ETr{(^)tf2*} < N 2t+2 a t+1 h+^^lt e ^j , (13) 
where C(a, k) is N- and t -independent. 



Let Cmax be the eigenvalue of f2 with the largest absolute value, and denote by 
•Sn(p) Q {+!> — 1}^^ the set of signature sequences, such that ICmaxI < NpX^\i°°\ 
Applying Markov inequality to the random variable |Cmax| 2 * < Tr it is 

easy to show that 
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By properly choosing p and t (and setting Sjy = Sn(p)), this in turns imply 
P {Sn} < 0(N- a ) with a > 0. 

We must now prove that, for any signature sequence in Sn, = [T^J — > 0. 
Since A(s) -» A(°°\ A(s) -» A(°°), there exists s*, such that |C max | < N p\(s*)X(s*), 
with p < 1. Defining Q = -Q/NX(s*)X(s*), we have 
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(15) 



Since the largest eigenvalue of has modulus smaller than one, 

Furthermore, since X(s), X(s) increases with s, the number in curly brackets is 

smaller than one. Therefore —* as well. 

Unlike for the means Xi, the BP estimate of variances is only correct in the 
large system limit N,K — > oo. In order to quantify the error made, we define 



(1/L 
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Theorem 2 Assume the S{ a 's are i.i.d. uniformly random in {+1, —1}, then D — > 
in probability as N — > oo. 

Sketch of proof: Freeman and Weiss [7] provide an explicit expression for the 
discrepancies <5j. In the present case, their expression can be shown to be equivalent 
to the following one 




Here the sum runs over all the closed non-reversing paths w, on the bipartite graph 
in Fig. ^ starting and ending at i. 

In order to convey the basic idea of the proof, let us restrict the above sum to 
the paths of a given fixed length 2t, and call 5® the corresponding quantity. It is 
easy to show that 

E ^) 2 = p4j F (^4^) 2 'lM|, d») 

where S l N K (2t) is the set of couples of closed paths starting at i on the bipartite 
graph in Fig.0 such that each edge in the graph is visited an even number of times. 
The dominating contribution to 1^/^(24)1 comes, in the large system limit, from 
couples of coincident paths. Their number scales like N 2t ~ l : we are free to choose 
each step but the last one among O(N) vertices. Therefore E(^*' ) ) 2 = 0(N~ X ) — > 0. 

In order to complete the proof, one has to control those paths in Eq. (|17|). whose 
length diverges with N. 



3 MMSE receiver recovered 

The per- iteration complexity of the BP algorithm considered above, scales N 3 . 
In fact the number of messages is equal to the number of edges in the complete 
bipartite graph of Fig. i.e. 0(N 2 ), and each message update involves O(N) 
sums. This may be too much for some applications. Also, although we proved that 
BP returns the correct MMSE estimate x, the relation between the two approaches 
is quite puzzling. Recall that MMSE is given by the explicit formula 



x 



(y) = [er J, I + S t S] -1 S t y. (19) 



where § is the TV by K matrix whose columns are the signature sequences of the 
users. This formula only involve vectors of size O(N). BP, on the other hand, 
involved linear operations on the vector py^J, of size NK. 

In order to clarify the relationship between the two approaches, let us write 



(t) ~G?+S4i a , £ = ^ +5 e. (20) 



The BP update equations for 7^ a , cf . Eq. © , read in the new variables (to lighten 
the formulae we assume here Si a G {+!,—!}) 



= (21) 



*r£? = -^-^s ia G^-^-s ia 5tL- (22) 
\{t)VN X(t)y/N 

Analogously, the update equations for become 

e °" = "■-a^S^-^S'*^- (23) 
5 ^ = mn"^ + mn s '^- (24) 

From Eqs. dHJ), it follows that S-ff^S-ff^ = 0{N- l l 2 ) and are given ap- 

proximately by 

*>" = -JfW"^" + 0(Ar,) ' ^ = ^vw""^ + 0(JV_1) ■ 

(25) 

Substituting in Eqs. (|2*Tj) . ((2*5|) . and neglecting 0(N~ 1 ) terms, we get, with a slight 
abuse of notation 



G? +l) = _4^Gf + YsaPP, (26) 



g(f) = y a + ^ Gt l) l -^Ys ka Gf. (27) 

A(t)A(t-l) Kt)^h 

These update equations will be referred to as 'approximate belief propagation'. As 
opposed to ordinary BP, they involve quantities associated to the vertices of the 
underlying graphical model. Their complexity per iteration is 0(N 2 ), which can 
be compared to the 0(A r3 ) complexity of ordinary BP for the present problem. 



Let us look for a fixed point G^ = [G^], G^ = [G^] of Eqs. ©, (J27J) 
and neglect (^(A^ 1 ) terms for the time being. Recall that, at the fixed point 
A(°°) = 1 + A, and A (oo) = 1/A, where A solves Eq. @. Substituting in Eqs. 1(251) . 
1(2 7(1 . and adopting matrix notation we obtain the fixed point conditions 

£(00) = G (oo) + A gtg(oo) (28) 

1 + A 

g(oo) = ^_ G (oo) + L_§ G (°°). ( 29 ) 

1+A 1+A v ; 

As explained in Sec. Q the conditional expectations estimates are given by x^ 00 - 1 = 
GH/L(»). It is not hard to realize that Z,(°°) = A^ + CKA^ 1 ) = l + A + O^- 1 ). 
By solving Eqs. (J2BJ) and (J2JJ) for G (oo) and using this formula, one finally gets 
aK°°) = where x(y) is given by Eq. (|19|) . 

The above discussion has important practical implications. One can replace 
the original BP equations, with Eqs. (|26[) . (|27|) . which have lower complexity. This 
provides an 'approximate BP' algorithm. A more careful examination of the above 
steps allows to prove the following statement. 

Proposition 1 Assume the Si a 's are i.i.d uniformly random in { + 1,-1}. Then, 
with high probability, approximate BP is convergent, G^' — ► G^°°\ — > L^°°\ 
and provides correct estimates of the conditional means, i.e. = G^'/L^ 00 ' = 
x(y). 



4 Numerical simulations and convergence rate 

In the previous Section we discussed the complexity per iteration of the BP detec- 
tion algorithm. We found that it can be as small as 0(N 2 ) if the 'approximate' 
update equations (|26[). ((27|) are used. Optimal (MMSE) detection is achieved after 
convergence to the BP fixed point. 

One may wonder how many iterations are required for the algorithm to get 
'reasonably close' to its fixed point. A first answer is provided by the proof of 
Theorem ^ Each BP iteration involves the multiplication of the vector of messages 
by the matrix Q, and the division by NX(t)X(t). Lemma |2J suggests in turn that 
the eigenvalues of ft, have absolute value not larger than Ny/a. Therefore, at each 
iteration, the distance between current vector of messages and the the fixed point 
messages is, roughly speaking, rescaled by a factor 

^^-^4, (30) 
X(t)X(t) 1 + A 

where we approximate the quantities on the left hand side by their value at the 
fixed point. An equivalent way of formulating this result consists in defining 

U(a,a)=- (log 1 , (31) 

where A is understood to be the solution of Eq. Q. After t* = t*(a,a) iterations, 
the distance from the fixed point is rescaled by a constant factor e~ 1 . Any precision 

5 can be achieved within a number of iterations t(5) ~ t*(cr, a) log(A/<5). where A 
is the distance between the initial and fixed point messages. 

It can be shown that, the BP estimate after the first iteration is equal to the 
one provided by the traditional matched filter receiver. As the number of iterations 
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Figure 2: Mean square error for BP based multiuser detection, as a function of the 
number of iterations. The asymptotic values coincide with the ones for MMSE detection. 
Continuous lines refer to ordinary BP, and dashed lines to approximate BP. Left: K = 
100 users and N = 200 chips (a = 0.5). Right: K = 100 users and N = 100 chips 
(a = I). 

increases, BP moves from the matched filter to the MMSE estimate. After about 
£*(cr, a) log(l/e) iterations, the distance from the MMSE estimate is reduced by a 
factor e with respect to a matched filter receiver. 

In order to confirm this analysis and obtain some concrete feeling of the ac- 
tual number of required BP iterations, we plot in Fig. |2] the results of numerical 
simulations of the BP based detection algorithm. Each couple of curves refers to 
the application of BP (continuous curves) or approximate BP (dashed curves) to a 
single signature/noise realization. Convergence to the MMSE receiver was achieved 
in each case. 

A small number of iterations (about 5) provides a large improvement over the 
matched filter receiver. The number of iterations required for convergence is always 
moderate and consistent with the estimate of Eq. (j31j) . For a = 0.5 (left frame), 
we obtain, for instance, t* 2.7, 2.4, 1.7 and 1.0 for (respectively) a = 0.1, 0.2, 
0.4 and 0.8. If a = 1 (right frame), we obtain ~ 10.0, 5.0, 2.5 and 1.3 for the 
same values of a. 

5 Conclusion and generalizations 

We can envision two types of generalizations of the present work. First, one can 
refine the model Q by considering, for instance, different power constraints for 
different users, or signature sequences with more general distributions than the one 
considered in this paper. While such generalizations can be technically cumber- 
some, they should possible along the same lines exposed here. 

The second direction consists in modifying the a priori distribution of the signals 
Xi, by considering, for instance, binary antipodal signals, x,- t S {+1,-1}. The 
proves of Theorems ^ and El are heavily based on the general results on belief 
propagation for Gaussian graphical models in [7] and [Hj. These results do not 
extend to symbols x^s with a general a priori distribution. However, we expect 
the thesis (convergence and correctness of belief propagation) to remain true in the 
large system limit. We plan to use the ideas of (Hj to prove this claim. 
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